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ABSTRACT 

We investigate the cosmological evolution of GRBs, using the total gamma-ray 
fluence as a measure of the burst strength. This involves an understanding of 
the distributions of the spectral parameters of GRBs as well as the total fluence 
distribution - both of which are subject to detector selection effects. We present 
new non-parametric statistical techniques to account for these effects, and use these 
methods to estimate the true distribution of the peak of the vF v spectrum - E p - 
from the raw distribution. The distributions are obtained from four channel data 
and therefore are rough estimates; hence, we emphasize the methods and present 
qualitative results. Given its spectral parameters, we then calculate the total fluence 
for each burst, and compute its cumulative and differential distributions. We use 
these distributions to estimate the cosmological rate evolution of GRBs, for three 
cosmological models. Our two main conclusions are the following: 1) Given our 
estimates of the spectral parameters, we find that there may exist a significant number 
of high E p bursts that are not detected by BATSE, 2) We find a GRB co-moving 
rate density quite different from that of other extragalactic objects; in particular, it is 
different from the recently determined star formation rate. 

Subject headings: gamma rays: bursts - methods: statistical - cosmology: 
observations - miscellaneous. 

1. Introduction 



Identification of several gamma-ray bursts (hereafter GRBs) with X-ray, optical and radio 
afterglows, as well as measurements of the redshifts of GRB 970508, GRB 981214, and GRB 
980703 are beginning to define the GRB paradigm. However, this information is not sufficient for 
detailed cosmological studies of their distribution and evolution. Until more information on the 
distances to GRBs are revealed, we must continue to rely on standard statistical studies (such as 
obtaining number count distributions of fluence or flux) to gain further insight on the nature of 
these events. 



However, information on the spatial distribution - although necessary - will not be sufficient 
for solving the puzzle of GRBs. In particular, we need spectral information to understand the 
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energy release and radiation processes at the burst. This paper deals with the determination of 
the spatial and spectral properties of GRBs. 

The information on the spatial distribution comes primarily from the investigation of source 
counts or the logiV-logS' relation. Assuming a luminosity function (LF) and a cosmological 
model, one can convert such counts to spatial (or redshift) distributions (Weinberg, 1972). 
The likelihood of the assumed luminosity function determines our confidence in the redshift 
distribution or cosmological evolution of these sources. The usual practice is to use the peak flux 
f p (v) as a measure of the burst strength. This flux is related to the peak luminosity L p (y) via 
Lp{y) = f p (v)VLd 2 L K(z), where cIl(z) is the luminosity distance, K(z) is the so-called K-correction 
due to redshift of the spectrum, and f2 is the solid angle into which the radiation is beamed 
(Q = 4-7T for isotropic emission). However, we have no knowledge of the burst peak luminosity 
function - say ^(L p (v)) - and it is unclear what assumptions can be made about this function. For 
example, in a fireball model, both the duration and the peak luminosities are expected to depend 
on the details of the fireball (e.g. the bulk Lorentz factor and amount of baryon loading), so that 
their distributions will be broad and perhaps comparable to the observed range of the peak fluxes. 
As a result, the luminosity function cannot be readily de-convolved from the spatial distribution, 
and we cannot get an accurate measure of the burst evolution or spatial distribution from the 
logA^-logS 1 diagram. 

On the other hand, in a neutron star (or black hole) merger model, the total energy released 
- represented by the gravitational potential energy of the merging process - is expected to be well 
defined and have a narrow distribution; this is also most likely true for the "hypernova" model. 
Hence, we would like to determine what radiative signature is a representative measure of this 
energy and the strength of the GRB. 

As the name "gamma ray burst" implies, the peak gamma-ray flux f p (or luminsity L p ) of a 
GRB far exceeds the peak fluxes at other energies. However, until recently, we did not know at 
what frequency the total radiated energy peaked. The recent X-ray, optical and radio observations 
of the afterglows - in particular the fact that these fluxes decline more rapidly than \jt (see e.g. 
Murakami et al., 1997, Galama et al., 1997) indicate that the energy fluence in the gamma-ray 
range is also higher than in any other band. For isotropic emission or a beaming angle independent 
of wavelength, this implies that the total gamma-ray radiated energy, £ tot , is the best measure 
of the total released energy. Therefore, it is plausible that the distribution of the gamma-ray 
radiated energy is narrow, so that the total gamma-ray fluence Ftot = ^tot{^ + z )/{^^l) provides 
the best indication of the strength of the burst for the logA^-logS analysis and other investigations 
(assuming VL has a narrow distribution). In particular, it is a better measure than the observed 
fluence F Q b s , which represents the fluence in some narrow energy range (i.e. the energy range in 
which the detector triggers). 

Note that the relation between the monochromatic flux (or fluence) and the corresponding 
luminosities require a knowledge of the so called K-correction due to redshift of the spectrum. 
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This complication is not present in the relation between the total radiated energy and the total 
gamma ray fluence. It is also important to note that F to t depends on the spectral properties of 
the burst. Hence, if we want to simplify matters by dealing with the total fluence, we need to 
understand the spectral distribution F(y) of the GRBs. 

The primary goal of this paper is to demonstrate how we can obtain information on the 
cosmological evolution of GRBs, using the observed distribution of F to t- However, for an accurate 
determination of the distribution of F tot , we need to determine the spectral parameters of each 
burst. It is common to use the following four parameters (say, in a Band spectrum (1993)) to 
characterize a burst's spectrum: the photon energy at the peak of the vF v spectrum, E p , the low 
energy spectral index, a, the high energy spectral index, (3, and a normalization A. The best 
source of data for the purpose of determining these parameters and F tot is the BATSE instrument 
on board the Compton Gamma- Ray Observatory (CGRO). However, because the BATSE trigger 
criterion is based on photon counts between 50 - 300 keV (for the 3B catalog), the selection biases 
or thresholds on F to t as well as the spectral parameters can be quite complicated. In §2 we show 
how the trigger criterion used by BATSE can be used to obtain the selection bias on any measured 
quantity - in particular, Ftot an d E p - and describe new methods to account for the resultant 
complicated data truncation. In §3, we present our results on the determination of the spectral 
parameters. The distribution of E p is given in §4, and in §5 we derive the distribution of F tat , give 
a scenario for cosmological evolutions of the GRBs, and discuss its relation to the new observations 
of afterglows and evolution of other cosmological sources. A brief summary is presented in §6. 

2. Description of Data and The Basic Problem 
2.1. The Data 

BATSE has produced a large catalog of GRBs with a wealth of information on their temporal 
and spectral characteristics. This detector has been designed to optimally detect GRBs and 
has been very successful in achieving this goal. Nevertheless, like all detectors, it has intrinsic 
limitations and may miss certain types of bursts or provide only limited information on others. 
BATSE is a triggered instrument; it will record data on a burst if the detector counts in a time 
interval At and energy range AE exceed some threshold value C m i n (At) determined by the 
background fluctuations. BATSE uses three trigger time intervals At = 1024, 256, and 64 ms with 
four energy channels AE = 25 - 50, 50 - 100, 100 - 300, and > 300 keV. For most of the data, 
triggering is based on counts in channels 2 and 3; recently, however, other channel combinations 
have been tried. Here we will use the 4 channel LAD fluence values published in the BATSE 
3B catalog (Meegan et al., 1995), which utilizes the former trigger scheme, so that the observed 
fluence F Q b s = ^QkeV^ F(v)du. However, the formalism described below is applicable to the latter 
data as well. Note that we analyze bursts with C max / C m i n on 1024ms timescale. Because of the 
lack of resolution of the four channel data, here we will emphasize the method of our calculations 
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rather than give concrete quantitative results. Future work will involve calculations made from 16 
channel CONT data. 



2.2. Data Truncation 



The BATSE detector trigger conditions lead to several limitations on the data. The most 
obvious is that only bursts with peak counts C max > C m i n are detected. Determination of the 
average V/V max = (C m j„/C mai ) 3 / 2 , and its distribution were one of the primary objectives of 
BATSE using this trigger configuration (Meegan et al., 1992). Another limitation arises from the 
finite durations At of the triggering criteria. As pointed out by Petrosian and Lee (1994, see also 
Lee and Petrosian, 1996), this reduces the detector efficiency of short duration bursts (duration 
T < At) in the sense that there is a relative bias against detection of weaker bursts with shorter 
duration. Similarly, as pointed out by Piran and Narayan (1996), the finite range AE of the 
detector energy channels provides another limitation (see also Cohen, Piran, and Narayan, 1998). 
In this case, the detector efficiency is highest for bursts with most of their photons in the range 
AE, and decreases for bursts with more of their flux outside this range. 

It turns out that we can account for these limitations by utilizing the BATSE trigger 
condition: Let us consider a burst characteristic, say X, whose value is measured by BATSE. For 
each burst, then, we know X, C max , and C m i n . Now, in the spirit of the V/V max test, one can ask 
what is the possible range of values of X that this burst can have and still trigger the instrument. 
In general, X can have an upper limit u, a lower limit /, or both limits; X G T = [l,u]. This 
question was first posed by Petrosian and Lee (1996) in connection with the determination of the 
detection threshold or efficiency of the observed energy fluence F f, s ; in this case, the fluence has 
only a lower limit F^sUm- It was shown that this threshold (or for that matter, the threshold on 
any other measure of burst strength which is proportional to the photon counts) is obtained from 
the simple relation 

Fobs I Fobs ,lim = C max /Cmini (1) 

provided the burst's spectrum does not change drastically throughout its duration. A similiar 
relation holds for the total fluence F tot as we describe in §5. 

For other bursts characteristics, this may not be the case. For example, we have shown (Lloyd 
and Petrosian 1998) that the spectral parameter E p has both an upper and lower limit. The 
values of E Pimax and E Ptm i n are not related to C max /C m i n and E p in a simple form as above, and 
thus we require a more complex procedure to determine their values. We will describe this in §3. 
Given the thresholds, we can then use non-parametric statistical techniques to obtain an estimate 
the parent distribution of the relevant observable from the observed distribution. 



2.3. Dealing with Data Truncation 
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Correcting for data truncation in the one-sided case was discussed in detail by Efron and 
Petrosian (1992). This method has been used to determine the correlation between parameters 
in a bivariate (or multivariate) setting, and to obtain univariate distributions of several GRB 
characteristics (see, e.g., Lee and Petrosian, 1996). The basic idea is that one can use information 
from the untruncated region of a distribution to estimate how much data is missed due to the 
truncation. For a description of the basic process, see also Petrosian (1992 or 1993). 

Dealing with data which has both upper and lower limits is more complex and requires a 
generalization of the above technique. Recently, Efron and Petrosian (1998) have developed a 
method to deal with this kind of data truncation. We briefly review this method here. 

Consider data points (xi,yi) (in our case (E p , F b s )), where Xi has lower and upper limits li 
and Ui, respectively; Xj G Tj = [k,Ui\. Let us also assume that Xi and yi are uncorrelated. Let 
g(x) be the true distribution of x, which would be observed if there were no truncations. However, 
because X{ is limited to Tj, we observe the conditional distribution g(xi\Ti) = g(xi)/G(Ti), where 
G(Ti) = [J2j9( x j) '■ x j £ Ti\ is the probabilitly that x exists in Tj. We normalize g(xi) such that 
J2iLi <?i = 1> an d define the vectors 



g i = g{x i ) and Gi = G{T$ (2) 

and the matrix 



Jij 



"10, xjtTi. {6) 

Furthermore, by definition, one can show that Gi = Jijgj (where we have used the convention that 
repeated indices are summed over). The goal is to estimate g(x) from k,Ui, and Xi assuming all N 
cases are independently distributed. It can be shown that maximizing the likelihood (see Maloney 
et al., 1998) of the observed data set gives 

l/Si = Jji * (l/Gj) + const, (4) 

where Jji is the transpose of Jij . The procedure amounts to solving this equation iteratively, under 
the normalization condition and definition of Gj and J^ given above. Convergence is reached 
when const in the above equation goes to zero. This method can be used to determine univariate 
cumulative (G) and differential (g) distributions. As described in Efron and Petrosian (1998), it 
can also be used to determine correlations between relevant variables (e.g. fluence and E p or peak 
flux and E p ). 

We will use the methods described above to correct for data truncations and produce bias 
free distributions of the peak energy E p and the total fluence F to t ', we will then directly use the 
latter distribution to investigate the cosmological evolution of GRBs. 



3. Spectral Characteristics 



-6- 



The spectral properties of GRBs provide the most direct information about the physical 
processes associated with the event. In particular, the distribution of the spectral parameters 
and the correlation between these and other GRB characteristics can shed significant light on the 
radiation mechanisms and energy production of the burst. 

Therefore, our first task is to determine the spectral parameters of each burst. We characterize 
the spectral forms of the bursts by four parameters - a and (3 (the photon power law index at 
low and high energies respectively), E a (a break energy related to the peak of the vF v spectrum) 
and A (a normalization in units of ergs /cm 2 ). We choose two functional forms to represent the 
energy fluence spectra - a Band spectrum, and a generic double power-law spectrum with a 
smooth transition between the low energy and high energy power law behavior: 

UA/E ){E/E o r+^[{p-a){E/E -l)\, E < E 
t{tj) \(A/E )(E/E Y+\ E>E [b) 

F(E) _ {A/E ){E/E )Y^ 

F{E) ~ 1 + {E/E )°-P (6) 
We are interested in bursts whose vF v spectra have a maximum, so that we may correctly speak 
of E p as the "peak of vF v ". To define E p (as well as to keep the burst's total fluence Ftot finite), 
we require a to be greater than -2 and (3 to be less than -2 (which also implies a > (3). As we 
discuss below, this does not impose any large bias in our sample. Then, with these constraints, we 
can define E p for each spectrum above: E p = E {^^) and E p = E ( ~ ( f^ )~^" f° r (5) an d (6), 
respectively. 

From four channel data, we can determine the four spectral parameters with no degrees of 
freedom. However, this introduces strong interplay between E p and (3 or a. We avoid this by 
selecting (3 from an assumed gaussian distribution with a mean of k, —3, and a standard deviation 
of 1 (and of course the constraint (3 < —2). This allows for one degree of freedom. In fact, we 
have performed this analysis using several distributions for (3 (i.e. a uniform distribution, or 
a gaussian with a different mean and standard deviation), and have found similiar qualitative 
results. However, we should note that steeper, more negative /3's will naturally produce higher 
E p s, to accommodate the fluence in high energy channels (3 and 4). We find the rest of the 
parameters by minimizing x 2 via a downhill simplex routine (e.g. see Press et al., 1992). For 
the Band spectrum, we find acceptable fits for 291 bursts, while 268 give acceptable fits for the 
second spectrum listed above, out of a total of 518 bursts. Figures 1(a) and 1(b) show the raw 
distributions of the spectral parameters obtained with this method for the two spectral forms in 
equation 5 and 6, respectively. 

It turns out that the distribution of the parameters for bursts with unacceptably high x 2 is 
not very different from those with acceptable x 2 ■ We have tested to see if the reason for many high 
X 2 bursts came from the constraints on a and (3 mentioned above. Relaxing these restrictions, we 
found that only ~ 1% of bursts that gave acceptable fits had an a < —2, (3 > —2, or a > (3. We 
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therefore assume that the bursts with good fits are representative of the total sample. However, 
because the lack of energy resolution of the four channel data, we would like to emphasize that 
our fits are rough estimates of the actual spectral parameters, so that the following analysis 
exhibits our method rather than presenting a final quantitative result. As mentioned previously, 
in future work we will use fits from 16 channel CONT data, which will provide a more accurate 
determination of these parameters. 

Several studies have obtained the spectral parameters for the brightest GRBs (e.g., Band et al. 
1993) and investigated correlations between these parameters and peak flux, duration and spatial 
distributions (Mallozzi et al. 1995, Kouvelioutou et al. 1995, Pendleton et al. 1996). However, 
caution is required in such studies, when dealing with the raw data. The effects of selection biases 
must be evaluated before obtaining these distributions and correlations. The parent distributions 
of the spectral parameters could be quite different from the observed distributions shown in 
the above figures. The differences are caused by the truncations resulting from the triggering 
procedure, and are a measure of trigger efficiencies for each of the four spectral parameters 
a,P,E p , and A or F to t = A£(a,(3). The truncation on A or F to t is described by equation (1): 
Ftot,Um = Ftot (Cmin/C max)- As mentioned previously and shown in detail in the next section, 
there are also selection biases against high and low values of E p . For instance, BATSE is most 
likely to see bursts with E p in the energy range in which the detector triggers. Becuase the 
detector needs some minimum number of counts Cu m to trigger, then for each burst one can find 
an interval [E Pmin , E Pmax ] to which E p is confined in order for the burst to have been detected. 

Similarly, truncation is present for a and (3. For example, if a burst with given E p and A had 
a value of a which was postitive and large or had a value of (3 which was negative and large, - 
i.e. if the vF v spectrum had a steep rise or fall - then the observed C max could fall below the 
threshold C m i n and render the burst undetectable. Using the methods described in §2 and below, 
we can correct the raw distribution of all four parameters and determine their correlations with 
other GRB characteristics. However, the effect of the truncation on a and (3 is more complicated 
and not as pronounced as that of E p and F to t- We will deal with the former in future publications. 
Here, we limit our discussion to obtaining the true distributions of E p and F to t- We will present 
the results from the Band spectrum fits to the bursts only; the smooth double power-law spectrum 
in equation (6) gives qualitatively similiar results for all of the following analysis. 

4. Distribution of E p 

For each burst we have the values of a, (3, A, E p , and the observed fluence 
F obs = J** F a ^ tA (E,E p )dE, with E x = 50/eeU and E 2 = 300/ceU. From equation 1, we 
can calculate the fluence threshold -F f, Si zi m . Using this limit, we can determine the possible range 
of values E p can take on, and still trigger the BATSE instrument. It is clear that the limiting 
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value E Pt n m must satisfy the equation 

F a ,i3,A(E,E Pt u m )dE. (7) 

For the spectral forms and parameter ranges discussed in the previous section, the right hand 
side of equation (7) increases monotonically with E p n m , reaches a maximum for E p n m within 
the limits of the integration, and then decreases monotonically. Consequently, there will be two 
values of E Pt n m that satisfy this equation. To find these two solutions, we start with the value 
of E p determined by the simplex routine (for which the right side is equal to F b s ), and increase 
it (while keeping a, /?, and A equal to their determined values) until equation (7) is satisfied - 
that is, until the observed fluence is brought below the threshold. This value is the upper limit 
to Ep, E p ^ max . We then decrease E p until this equation is again satisfied. This gives the lower 
limit Ep yTnin . In essence, the procedure amounts to redshifting/blueshifting the burst until it is 
undetectable by BATSE. Figure 2 shows E p vs. E Pjmax and E Pt7n i n for the Band spectrum. Notice 
the truncation is much more prominent on the upper than the lower end; hence, we expect to see 
a more significant correction to the high end of the E p distribution. 

Using the method described in §2.3 for doubly truncated data, we corrected the observed 
distribution given that the observations can detect only bursts with E p limited to the interval 
[E Pmin , E Pmax ]. Figure 3 shows the raw and corrected differential distribution of E p for the 
Band spectrum. The raw and corrected distributions are normalized such that they are equal 
at low values of E p (where the data truncation is inconsequential). This normalization visually 
emphasizes the differences of the distributions at high E p . The figures suggest that a large number 
of high E p bursts are missed by the triggering procedure. 

In order to determine the significance of the difference between the raw and corrected 
distributions, we carry out the following two checks: First, we perform a K-S test on the 
cumulative distributions of E p , and a \ 2 test on their differential distributions for both the Band 
spectrum and the smooth double power law. For the \ 2 test, we divide the distributions into two 
bins, with equal numbers of observed bursts per bin. As shown in the table below, in both cases 
the tests indicate the probability that the observed and corrected E p distributions are the same is 
extremely small. 



Test 


Spectrum 


Probability the distributions same 


K-S 


Band 


1.8 x 1(T 6 


x 2 


Band 


3.0 x 1(T 7 


K-S 


Double Power Law 


4.4 x 1(T 5 


x 2 


Double Power Law 


2.4 x 1(T 5 



Second, we have carried out simulations to determine the accuracy of the method described in 
2.3; Appendix A contains these results. We begin with a parent distribution of E p 's, simulate an 
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observational truncation (defined by the limiting fluence of the burst; see §4). From this we obtain 
our observed distribution of E p 7 s. We then apply the above technique to correct for the truncation, 
and retrieve our parent distribution to very high accuracy, demonstrating the robustness of the 
method (see Appendix A for more details). 

Our results are in qualitative agreement with the observational data from the Solar Maximum 
Mission (SMM) presented by Harris and Share (1998a, 1998b); SMM is sensitive to higher energies 
than BATSE, so that it suffers the bias against detection of bursts with hard spectra to a much 
lesser degree. The results presented in Harris and Share (1998b) agree qualitatively with our 
corrected E p distribution. We perform a similiar x 2 as above, comparing the E p distribution from 
SMM with both the uncorrected and corrected BATSE E p distributions presented in Figure 3. 
In this case, we divide the data into two bins below and above 550 keV. We find the probability 
that the E p distribution from SMM data is the same as the raw BATSE E p distribution to be 
1.3 x 1CT 6 . On the other hand, we find a 90% probability that the E p distribution from the SMM 
data is the same as the corrected BATSE E p distribution. Note that the statistical method we 
use can correct the E p distribution only for the values of E p observed by BATSE. Therefore, we 
cannot determine the distribution above ~ 1 MeV (where the raw BATSE distribution ends). 
Similarly, because we have no observed E p s below ~ 20 keV, we cannot confirm Strohmeyer et al. 
(1998) reporting a significant number of low E p bursts (< 20 keV) observed by GINGA, but not 
observed by BATSE. 

A bias against detection of high E p bursts has significant implications on correlations of E p 
with other measures of burst strength like peak flux, F b. s , or F tot . The detector is most likely to 
miss the high E p bursts with low strength, so that a simple correlation analysis between raw values 
of E p and burst strength without the consideration of the data truncation can lead to erroneous 
correlations. Examples of studies which may be affected by such a bias include the correlation 
studies of Nemiroff et al. (1994) and Mallozzi et al. (1995). This will also have an important 
effect on determining whether or not the correlations are intrinsic or due to cosmological redshift 
(Brainerd, 1997). The above results indicate that a more thorough analysis of the correlations is 
required before conclusions can be reached on the redshifts or intrinsic effects. We will address this 
question in a future publication, when we have more accurate estimates of the spectral parameters. 

5. Distribution of the Total Gamma Ray Fluence 

As discussed in the introduction, an important way to learn about the spatial distribution 
of GRBs is to study the distribution of some standard candle variable of the event. We also 
conjecture the total radiated energy 8 tot in the gamma-ray range may be the best candidate for 
such a standard candle variable, and that the total fluence Ff t of the GRB may provide the best 
measure of distance. Consequently, we study the distribution of the bursts' total fluence; if the 
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burst radiates isotropically, the total fluence is related to the total radiated energy £ tot as: 



£tot 



(1 + *), (8) 



where d,L, the luminosity distance, depends on redshift and the cosmological model (Weinberg, 
1972). 

Again, because the detector is sensitive only over a finite energy range, BATSE does not 
necessarily obtain all photons from the burst. Hence, the detector measures only a portion 
of the burst's total fluence. However, the total fluence of a GRB can be obtained from the 
spectral fits to the observed fluence simply by integrating over the spectrum of that burst: 
Ftot = J*£^ in F aj p : A(E, E p )dE., where E m - m denotes the beginning of the gamma ray energy range 
(as an estimate, we use E min = 25 keV). This fluence will have an observational lower limit in the 
same way that the fluence from 50 — 300 keV has, and its threshold is given by the generalization 
of equation (1): F tot j im = ^, min F tot . Following the steps described in Petrosian and Lee (1996), we 
first test for any correlation between Ftot and F tot n m and parametrically remove this correlation. 
We then use the method described in Efron and Petrosian (1994) for one sided truncated data (this 
is equivalent to the method of §2.2, taking the upper limit u to oo), and obtain the cumulative 
and differential distributions for the total fluence. 

Figures 4(a) and 4(b) show the cumulative and differential distributions of the total and 
observed fluence for our sample. Marked on each curve are GRB 970508 (circle), GRB 971214 
(cross), and GRB 980703 (triangle). Note that we do not find a good spectral fit for GRB 980703; 
hence, we obtain a lower limit to its total fluence simply by summing up the published fluence 
values for the four channels of the LAD detector. Following the tradition of radio astronomers (see 
Lee and Petrosian, 1996), we divide out the —3/2 power law dependence of the number on the total 
fluence. In this way, deviation from the horizontal line suggests deviation from a homogeneous, 
isotropic, static, Euclidean distribution. In these figures, we have included all 518 bursts available 
- including those whose spectral fits gave high values of x 2 - so as n °t to underestimate any part 
of these distributions. We have also applied the same procedure to the subsample of bursts with 
acceptable values of \ 2 an d find that the distributions are similiar to those of the complete sample 
within about 30 %. Note that the fluence distributions in Figures 4(a) and 4(b) are unlike the 
counts of other well known extragalactic sources such as galaxies, radio sources and AGNs or 
quasars. The transition from 3/2 power law is too abrupt and the slope beyond this transition is 
nearly constant. Clearly, some extraordinary evolutionary processes are at work. We explore this 
in more detail in the next section. 



5.1. Rate Evolution 



There have been several detailed analyses of the logA^-logS distribution for the peak fluxes of 
GRBs with inconclusive results. This is primarily because any observed distribution can be fitted 
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to an arbitrary luminosity function and evolution even if one assumes a cosmological model (see 
e.g. Rutledge et al. 1995, Reichert and Meszaros, 1997, Meszaros and Meszaros, 1995, Hakkila et 
al., 1996). Additional uncertainty associated with these results comes from neglect of the time 
bias, spectral bias, and uncertainty in the value of the K-correction. These source of error are 
absent when dealing with the total fluence, which requires no correction for time bias (Lee and 
Petrosian, 1996), includes the spectral bias (described above), and requires no K-correction. Thus, 
the relation between the total fluence counts and either the cosmological models or the luminosity 
function and its evolution is more straight forward. Here we concentrate on the rate evolution of 
GRBs. 

To obtain some indication of possible evolution, we assume several representative values for 
the total radiated energy £ tot and derive the comoving rate density from the observed differential 
counts of the fluences as 

p{z) = (1 + z)n(F tot )(dV/dz)-\dF tot /dz) (9) 

where n(F tot ) is the differential distribution of the total fluence. Note that in all of our calculations, 
we have assumed H Q = 60km s _1 Mpc _1 

Figures 5(a), 5(b) and 5(c) show p(z) for three representative cosmological models. In order 
to span the range of possible evolutions, we choose three extreme models: a flat matter dominated 
model (O m = 1, = 0), a curvature dominated model (Q m = 0.2, 17a = 0), and a flat vacuum 
energy dominated model (fi m = 0, = 1). Here, J7 m is the density parameter (equal to the 
ratio of matter density to the critical density), and 17 a = Ac 2 /3H% (where A is the cosmological 
constant). Within each figure, p(z) is plotted for four standard candle energies: Etot = 10 51 , 10 52 , 
10 53 , 10 54 ergs. Superposed on each figure is the star formation rate (SFR) from Madau et. al 
(1997, solid line), this rate delayed by 2 x 10 9 years (dashed-dot line), as well as the SFR convolved 
with a distribution of time delays P(t) oc t^ 1 (dotted curve), which may arise in a neutron star or 
black hole merger scenario (see e.g. Totani, 1997). Note that the SFR was calculated using an 
Q m = 1, $7 a = cosmology (Madau, 1997); however, the shape of this curve as well as the 2 other 
curves derived from the SFR are fairly insensitive to the cosmological model (see e.g., Totani, 
1997). 

Also marked on each figure are the rate densities of GRB 970508 (circle), GRB 971214 
(cross) and GRB 980703 (triangle), given their measured redshifts, and the energy £ tot required to 
produce a burst at this redshift (given its total fluence). Table II lists this required energy. Again, 
because we have only obtained a lower limit to the total fluence of GRB 980703, we only list a 
lower limit to its required radiated energy. 
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In all cases, the curve corresponding to Stot = 10 ergs is ruled out by the observational 
data. The £ to t > 10 52 ergs curve accommodates GRB 970508, while GRB 971214 and GRB 
980703 require £ tot > 10 53 ergs. In the Q A = 1 case, the required £ to t > 10 54 ergs for GRB 971214 
begins to exceed theoretical limits, which may be evidence against this cosmological model or 
evidence for strong beaming of the burst radiation. Our results suggest that the total gamma-ray 
radiated energy is not a standard candle and is spread over at least one decade. We emphasize 
that these results depend critically on the total fluence, which in turn depends on accurate values 
for the spectral parameters for these bursts. Higher resolution fits to the data as well as more 
measurements of the redshifts of bursts will allow us to constrain the curves more definitely. 

Finally, note that none of the three curves for p(z) follow the star formation rate curve, as 
one would expect in many GRB models such as the merger or hypernovae models. The density 
evolution for GRBS in unlike that for AGN's or ordinary galaxies as well. The results agree - at 
least qualitatively - with those of Totani (1998). However, this is in contradiction with recent 
claims that the GRB rate follows the star formation rate (Wijers et al., 1998). 

6. Summary and Conclusion 

The primary aim of this paper has been to investigate the cosmological evolution of GRBs, 
which involves understanding the distributions of the spectral parameters as well as the total 
fluence of the burst. Selection effects limit the information we can observe to explore these 
distributions. 

We have presented new non-parametric methods to account for these effects. We have applied 
these methods first to the distribution of the peak of the vF v spectrum, E p . We find that this 
spectral characteristic suffers both upper and lower thresholds, although the upper threshold has 
the dominant effect. We correct the observed distribution for this truncation and present the 
parent distribution of E p . This contains a large number of high E p bursts not evident in the 
observed distribution, which is in qualitative agreement with the SMM results of Harris and Share 
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(1998b). This is important for studying correlations between E p and other burst characteristics, 
as well as for theoretical considerations concerning GRB models. 

Using the spectral fits for each burst, we determine a rough measure of the total gamma-ray 
fluence and present the GRB source counts based on this measure. We obtain both the differential 
and cumulative counts of the total and observed fluences. 

We convert these distributions to the comoving rate density of GRBs for 3 different 
cosmological models and 4 different values of the total gamma ray radiated energies. We find 
that none of the curves follow the star formation rate as predicted by various GRB models. More 
observations on the redshifts of GRBs as well as high resolution spectral data (to obtain a precise 
value for the total fluence) are needed to substantiate these results. 

As a final note, we summarize some of the caveats in our analysis. First, the spectral 
parameters are determined based on four channel data and are only rough estimates of the true 
spectral properties of the burst. In particular, because (3 is chosen from an assumed distribution, 
the correction to the E p distribution could be overemphasized. These spectral parameters also 
affect the values of the total fluence (since Ftot = J F a p A{E, E p )dE), which in turn affect the 
density evolution functions for the bursts. Furthermore, to derive the density curves, we have 
assumed that the bursts total fluence is a standard candle; although this seems to be the most 
plausible candidate for a standard candle, the observations of GRB 970508, 971214, and GRB 
980703 indicate the contrary. Note, however, that we can use a standard candle upper and lower 
limit to at least constrain the possible evolution of GRBs. Finally, effects such as beaming were 
not taken into account in the density evolution analysis. For a burst beamed into a solid angle £1, 
the total radiated energy 8 tot will decrease by a factor of Q/An from the isotropic total radiated 
energy; hence, the observed dispersion in the total radiated energy (about an order of magnitude 
or higher) could be due to variation in the degree of beaming. Meanwhile, beaming causes the 
number of burst occurances over some time interval to increase by a factor of 47r/f2. We will 
address the issue of how beaming can affect the cosmological rate evolution in a future publication. 

We would like to acknowledge P. Madau for the star formation rate data used to generate 
the three SFR curves in 5(a), 5(b), and 5(c). We would also like to thank B. Efron for help in 
implementing the statistical methods discussed in the text. Finally, N.M. Lloyd would like to 
thank members of the BATSE team for many useful discussions, and their hospitality during her 
visit to MSFC. 

A. Appendix 

Here, we describe the results of simulations which exhibit the technique and demonstrate the 
accuracy of the procudure used for correcting doubly truncated data. For this analysis, we use 
the Band spectrum displayed in equation 5. We pull E p from a uniform distribution between 5 
and 35 keV. a is also drawn from a uniform distribution ranging between -0.9 and 5.9, as was f) 
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which ranged from -9.0 to -2.1. Finally, our normalization parameter A is taken from a uniform 
distribution from 12 x 10~ 9 to 18 x 10 _9 ergs/cm 2 . Note that we are not trying to simulate the 
real distribution of any of these parameters - we choose these distributions for the purpose of 
displaying selection effects and to what degree our method can account for them. 

Once we have the spectral parameters, we can calculate the fluence in each of the four BATSE 
LAD energy channels. The limiting fluence of each burst is chosen from a gaussian distribution 
with a mean of 6 x lCr n ergs/cm 2 and a standard deviation of 6 x 10 _10 ergs/cm 2 (note again, this 
is not an attempt to simulate the actual BATSE detector response). We then ask the question: 
How many of the 1000 simulated bursts have their fluence in channels 2 and 3, greater than the 
limiting fluence? We find 554 bursts survived this cut - we will call this the observed sample. 
We then proceed to apply the method described in §2 to our observed sample, to see if we could 
get back the original distribution of E^s. Figure 6 shows the the parent sample, the observed 
sample, and the correction to the observed sample. As evident, the correction accounts for the 
truncated data, and approximately reproduces the original sample. See §2 for a full description of 
the technique. 
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Fig. 1. — The distribution of the spectral parameters determined from a) the Band spectrum and b) 
a smooth double power-law spectrum respectively. Note that (3 is drawn from the same truncated 
gaussian distribution in both cases, and is not a fitted parameter. 




Fig. 2. — The maximum and minimum values of E p for the Band spectrum. Note that the 
truncation is more severe from above than below, which indicates that observations may miss 
a population of GRBs with high E p . 
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Fig. 3. — The observed and "corrected" distribution of E p for the Band spectrum. The figure 
suggests that there is a significant sample of high E p bursts undetected by BATSE. As mentioned 
in the text, the methods we used will not produce a correction beyond the raw observed distribution. 




Fig. 4. — The cumulative (a) and differential (b) distribution of both the observed (F b s ) and total 
(F tot ) fluences. The circle indicates the observed and total fluences of GRB 970508, the cross marks 
GRB 971214, and the triangle marks GRB 980703. In each figure, note the similarities between 
the two distributions; they both display an abrupt break from the HISE dependence (see text), 
indicating presence of strong evolutionary processes. 
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Fig. 5. — The comoving rate density versus redshift z for 3 different cosmological models indicated 
by the values of f2 m and Q\. Within each plot, we show the rate density for 4 different standard 
candle energies; the numbers 0.1, 1, 10, and 100 mark the total gamma-ray radiated energy in 
units of 10 52 ergs. Note that these curves are significantly different from: 1) the star formation rate 
(SFR) from Madau et al., 1997 (solid line), 2) the star formation rate with a time delay of 2 x 10 9 yr 
(dot-dash line) , and 3) the star formation rate convolved with a distribution of time delays (dotted 
line, Totani et al., 1997). GRB 970508 (circle), GRB 971214 (cross) and GRB 980703 (triangle) 
are indicated, given their measured redshifts and the energy Etot required to produce a burst at 
this redshift (given its total fluence, see text). Note that no single standard candle energy can 
accomadate the observed redshifts of these bursts, in any of the three models. 
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Fig. 6. — The original E p distribution (solid line), the observed E p distribution given the criterion 
Fobs > Fiim (dashed line), and the corrected distribution using the method described in §2 of the text 
(dotted line). These distributions are obtained from the simulations described in the Appendix. 
Note that the procedure recovers the original distribution almost exactly. Any differences are 
primarily due to differences in the binnings of the distributions. 



